/**
 * 这个文件定义了基于sequence的一些操作
 * 包括: filter, mapReduce, pack, pack2, scan, split, sum等等
 * 其中后缀名带serial的都是串行实现，不应该被外部调用。
 * 不带serial后缀的方法可以被外部调用,这一类方法都做了granularity control.（较小的数据会调用串行算法）
 *  * 
 * ********这个文件内容值得弄懂********
 * 1. filter 比较简单，depth是O(logn)。
 * 2. a. pack: pack takes an input a range a boolean indicator sequence, 
 * and returns a new sequence consisting of the elements of the range 
 * such that the element in the corresponding position in the indicator sequence is true.
 *    b. pack2:将输入依照两个indicator sequence打包为两部分，如果indicator sequence是一个boolvector，则分为false和true打包
 * 3.mapReduce： 先讲每个数据map，然后reduce
 * 4，scan：计算前缀和（注意区分prefixsun，在scan后，out[i]=in[0]+...+in[i])
 * 5. prefixSum: 计算前缀和（注意此处prefixsum实际上是prescan, out[i]=in[0]+...+in[i-1], out[0]=0, 返回为所有数的和
 */
#ifndef A_SEQUENCE_INCLUDED
#define A_SEQUENCE_INCLUDED
//#pragma once

#include <iostream>
#include "config.h"
#include "utils.h"

//使用了SIMD指令来快速计算bit位
// For fast popcount
#include <immintrin.h> //AVX .for _popcnt64
//#include <x86intrin.h>
#include <xmmintrin.h> //__SSE__
#include <mmintrin.h>  //__MMX__

using namespace std;

#define _BSIZE 2048
#define _SCAN_LOG_BSIZE 10
#define _SCAN_BSIZE (1 << _SCAN_LOG_BSIZE) //串行数据块的大小，即block_for的最小大小

template <class T>
struct _seq
{
  T *A;
  long n;
  _seq()
  {
    A = NULL;
    n = 0;
  }
  _seq(T *_A, long _n) : A(_A), n(_n) {}
  void del() { free(A); }
};

template <class E>
void brokenCompiler__(intT n, E *x, E v)
{
  parallel_for(intT i = 0; i < n; i++) x[i] = v;
}

template <class E>
static E *newArray(intT n, E v)
{
  E *x = (E *)malloc(n * sizeof(E));
  brokenCompiler__(n, x, v);
  return x;
}

//名字空间：sequence
namespace sequence
{

  template <class intT>
  struct boolGetA
  {
    bool *A;
    boolGetA(bool *AA) : A(AA) {}
    intT operator()(intT i) { return (intT)A[i]; }
  };

  //getA是一个函数对象，getA<ET,intT>(AA)(i) 访问的AA的第i个元素
  template <class ET, class intT>
  struct getA
  {
    ET *A;
    getA(ET *AA) : A(AA) {}
    ET operator()(intT i) { return A[i]; }
  };

  template <class IT, class OT, class intT, class F>
  struct getAF
  {
    IT *A;
    F f;
    getAF(IT *AA, F ff) : A(AA), f(ff) {}
    OT operator()(intT i) { return f(A[i]); }
  };

/**
 * nblock计算_n可以分为多少个大小为_bsize的block
 */
#define nblocks(_n, _bsize) (1 + ((_n)-1) / (_bsize))

/***
 * 大于阈值使用并行，否则用串行
 * 循环处理_st->ne之间的内容，如果_ne-st < thresh,使用并行循环
 * 
 * 注意：该方式似乎未被使用，更多的是用的block_for
 */
#define granular_for(_i, _st, _ne, _thresh, _body)  \
  {                                                 \
    if ((_ne - _st) > _thresh)                      \
    {                                               \
      {                                             \
        parallel_for(intT _i = _st; _i < _ne; _i++) \
        {                                           \
          _body                                     \
        }                                           \
      }                                             \
    }                                               \
    else                                            \
    {                                               \
      {                                             \
        for (intT _i = _st; _i < _ne; _i++)         \
        {                                           \
          _body                                     \
        }                                           \
      }                                             \
    }                                               \
  }
  /** 每次循环操作_bsize的内容
   * _i block的编号
 * _s block起始index ; _e block终止的index(不包含)
 * _bsize blocksize
 * _body 循环操作部分
 */
#define blocked_for(_i, _s, _e, _bsize, _body) \
  {                                            \
    intT _ss = _s;                             \
    intT _ee = _e;                             \
    intT _n = _ee - _ss;                       \
    intT _l = nblocks(_n, _bsize);             \
    parallel_for(intT _i = 0; _i < _l; _i++)   \
    {                                          \
      intT _s = _ss + _i * (_bsize);           \
      intT _e = min(_s + (_bsize), _ee);       \
      _body                                    \
    }                                          \
  }

  template <class OT, class intT, class F, class G>
  OT reduceSerial(intT s, intT e, F f, G g)
  {
    OT r = g(s);
    for (intT j = s + 1; j < e; j++)
      r = f(r, g(j));
    return r;
  }

  /**
   *  对i in [s,e)做reduce操作,reduce的算符g(a,b) 
   */
  template <class OT, class intT, class F, class G>
  OT reduce(intT s, intT e, F f, G g)
  {
    intT l = nblocks(e - s, _SCAN_BSIZE);
    if (l <= 1)
      return reduceSerial<OT>(s, e, f, g);
    OT *Sums = newA(OT, l);
    blocked_for(i, s, e, _SCAN_BSIZE,
                Sums[i] = reduceSerial<OT>(s, e, f, g););
    OT r = reduce<OT>((intT)0, l, f, getA<OT, intT>(Sums));
    free(Sums);
    return r;
  }

  template <class OT, class intT, class F>
  OT reduce(OT *A, intT n, F f)
  {
    return reduce<OT>((intT)0, n, f, getA<OT, intT>(A));
  }

  template <class OT, class intT, class F>
  OT reduce(OT *A, intT s, intT e, F f)
  {
    return reduce<OT>(s, e, f, getA<OT, intT>(A));
  }

  template <class OT, class intT>
  OT plusReduce(OT *A, intT n)
  {
    return reduce<OT>((intT)0, n, utils::addF<OT>(), getA<OT, intT>(A));
  }

  template <class intT>
  intT sum(bool *In, intT n)
  {
    return reduce<intT>((intT)0, n, utils::addF<intT>(), boolGetA<intT>(In));
  }

  //Map: g: IT->OT ,即g(A[i])是map函数
  //reduce: f: OT*OT->OT. 即f是一个reduce的二元函数，f(a,b)->r
  // g is the map function (applied to each element)
  // f is the reduce function
  // need to specify OT since it is not an argument
  template <class OT, class IT, class intT, class F, class G>
  OT mapReduce(IT *A, intT n, F f, G g)
  {
    return reduce<OT>((intT)0, n, f, getAF<IT, OT, intT, G>(A, g));
  }

  template <class ET, class intT, class F, class G>
  intT maxIndexSerial(intT s, intT e, F f, G g)
  {
    ET r = g(s);
    intT k = 0;
    for (intT j = s + 1; j < e; j++)
    {
      ET v = g(j);
      if (f(v, r))
      {
        r = v;
        k = j;
      }
    }
    return k;
  }

  /**
   * 
   * 
   */
  template <class ET, class intT, class F, class G>
  intT maxIndex(intT s, intT e, F f, G g)
  {
    intT l = nblocks(e - s, _SCAN_BSIZE);
    if (l <= 2)
      return maxIndexSerial<ET>(s, e, f, g);
    else
    {
      intT *Idx = newA(intT, l);
      blocked_for(i, s, e, _SCAN_BSIZE,
                  Idx[i] = maxIndexSerial<ET>(s, e, f, g););
      intT k = Idx[0];
      for (intT j = 1; j < l; j++)
        if (f(g(Idx[j]), g(k)))
          k = Idx[j];
      free(Idx);
      return k;
    }
  }

  template <class ET, class intT, class F>
  intT maxIndex(ET *A, intT n, F f)
  {
    return maxIndex<ET>((intT)0, n, f, getA<ET, intT>(A));
  }

  template <class ET, class intT, class F, class G>
  ET scanSerial(ET *Out, intT s, intT e, F f, G g, ET zero, bool inclusive, bool back)
  {
    ET r = zero;

    if (inclusive)
    {
      if (back)
        for (long i = e - 1; i >= s; i--)
          Out[i] = r = f(r, g(i));
      else
        for (intT i = s; i < e; i++)
          Out[i] = r = f(r, g(i));
    }
    else
    {
      if (back)
        for (long i = e - 1; i >= s; i--)
        {
          ET t = g(i);
          Out[i] = r;
          r = f(r, t);
        }
      else
        for (intT i = s; i < e; i++)
        {
          ET t = g(i);
          Out[i] = r;
          r = f(r, t);
        }
    }
    return r;
  }

  template <class ET, class intT, class F>
  ET scanSerial(ET *In, ET *Out, intT n, F f, ET zero)
  {
    return scanSerial(Out, (intT)0, n, f, getA<ET, intT>(In), zero, false, false);
  }

  // back indicates it runs in reverse direction
  template <class ET, class intT, class F, class G>
  ET scan(ET *Out, intT s, intT e, F f, G g, ET zero, bool inclusive, bool back)
  {
    intT n = e - s;
    intT l = nblocks(n, _SCAN_BSIZE);
    if (l <= 2)
      return scanSerial(Out, s, e, f, g, zero, inclusive, back);
    ET *Sums = newA(ET, nblocks(n, _SCAN_BSIZE));
    blocked_for(i, s, e, _SCAN_BSIZE,
                Sums[i] = reduceSerial<ET>(s, e, f, g););
    ET total = scan(Sums, (intT)0, l, f, getA<ET, intT>(Sums), zero, false, back);
    blocked_for(i, s, e, _SCAN_BSIZE,
                scanSerial(Out, s, e, f, g, Sums[i], inclusive, back););
    free(Sums);
    return total;
  }

  template <class ET, class intT, class F>
  ET scan(ET *In, ET *Out, intT n, F f, ET zero)
  {
    return scan(Out, (intT)0, n, f, getA<ET, intT>(In), zero, false, false);
  }

  template <class ET, class intT, class F>
  ET scanBack(ET *In, ET *Out, intT n, F f, ET zero)
  {
    return scan(Out, (intT)0, n, f, getA<ET, intT>(In), zero, false, true);
  }

  template <class ET, class intT, class F>
  ET scanI(ET *In, ET *Out, intT n, F f, ET zero)
  {
    return scan(Out, (intT)0, n, f, getA<ET, intT>(In), zero, true, false);
  }

  template <class ET, class intT, class F>
  ET scanIBack(ET *In, ET *Out, intT n, F f, ET zero)
  {
    return scan(Out, (intT)0, n, f, getA<ET, intT>(In), zero, true, true);
  }

  /*scan 操作,合并算子是plus，即out[i]=sum_{j=0}{i}In[i]*/
  template <class ET, class intT>
  ET plusScan(ET *In, ET *Out, intT n)
  {
    return scan(Out, (intT)0, n, utils::addF<ET>(), getA<ET, intT>(In),
                (ET)0, false, false);
  }

  template <class intT>
  intT enumerate(bool *In, intT *Out, intT n)
  {
    return scan(Out, (intT)0, n, utils::addF<intT>(), boolGetA<intT>(In),
                (intT)0, false, false);
  }

#define _F_BSIZE (2 * _SCAN_BSIZE) //blocksize in for iteration, _SCAN_BSIZE=2^10

  //返回bool数组F1的元素总和，F1的长度为n
  // sums a sequence of n boolean flags
  // an optimized version that sums blocks of 4 booleans by treating
  // them as an integer
  // Only optimized when n is a multiple of 512 and Fl is 4byte aligned
  template <class intT>
  intT sumFlagsSerial(bool *Fl, intT n)
  {
    intT r = 0;
    if (n >= 128 && (n & 511) == 0 && ((long)Fl & 3) == 0)
    {
      int *IFl = (int *)Fl;
      for (int k = 0; k < (n >> 9); k++)
      {
        int rr = 0;
        for (int j = 0; j < 128; j++)
          rr += IFl[j];
        r += (rr & 255) + ((rr >> 8) & 255) + ((rr >> 16) & 255) + ((rr >> 24) & 255);
        IFl += 128;
      }
    }
    else
      for (intT j = 0; j < n; j++)
        r += Fl[j];
    return r;
  }

  template <class intT>
  inline bool checkBit(long *Fl, intT i)
  {
    return Fl[i / 64] & ((long)1 << (i % 64));
  }

  template <class intT>
  intT sumBitFlagsSerial(long *Fl, intT s, intT e)
  {
    intT res = 0;
    while (s % 64 && s < e)
    {
      if (checkBit(Fl, s))
        ++res;
      s++;
    }
    if (s == e)
      return res;
    while (e % 64)
    {
      if (checkBit(Fl, e - 1))
        ++res;
      e--;
    }
    // Do the rest with popcount
    intT be = e / 64;
    intT bs = s / 64;
    for (intT i = bs; i < be; ++i)
    {
      res += _popcnt64(Fl[i]);
    }
    return res;
  }

  /**
   * 顺次打包f函数s到e的元素。
   * 即对于第i=[s,e),如果f(i)=true，将f(i)打包到out，否则忽略
  */
  template <class ET, class intT, class F>
  _seq<ET> packSerial(ET *Out, bool *Fl, intT s, intT e, F f)
  {
    if (Out == NULL)
    {
      intT m = sumFlagsSerial(Fl + s, e - s);
      Out = newA(ET, m);
    }
    intT k = 0;
    for (intT i = s; i < e; i++)
      if (Fl[i])
        Out[k++] = f(i);
    return _seq<ET>(Out, k);
  }

  template <class ET, class intT, class F>
  void packSerial01(ET *Out0, ET *Out1, long *Fl, intT s, intT e, F f)
  {
    if (Out0 == NULL)
    {
      intT m = e - s - sumBitFlagsSerial(Fl, e, s);
      Out0 = newA(ET, m);
    }
    if (Out1 == NULL)
    {
      intT m = sumBitFlagsSerial(Fl, e, s);
      Out1 = newA(ET, m);
    }
    intT k0 = 0;
    intT k1 = 0;
    for (intT i = s; i < e; ++i)
    {
      if (checkBit(Fl, i))
        Out1[k1++] = f(i);
      else
        Out0[k0++] = f(i);
    }
  }

  template <class ET, class intT, class F>
  void packSerial0(ET *Out, long *Fl, intT s, intT e, F f)
  {
    if (Out == NULL)
    {
      intT m = e - s - sumBitFlagsSerial(Fl, e, s);
      Out = newA(ET, m);
    }
    intT k = 0;
    for (intT i = s; i < e; i++)
    {
      if (!checkBit(Fl, i))
      {
        Out[k++] = f(i);
      }
    }
  }

  template <class ET, class intT, class F>
  void packSerial1(ET *Out, long *Fl, intT s, intT e, F f)
  {
    if (Out == NULL)
    {
      intT m = sumBitFlagsSerial(Fl, e, s);
      Out = newA(ET, m);
    }
    intT k = 0;
    for (intT i = s; i < e; i++)
      if (checkBit(Fl, i))
        Out[k++] = f(i);
  }

  /**
   *  计算前缀和，注意第i项的和不包括第i项。初始项为0
   *  即data[i]=data[0]+...data[i-1]
   */
  template <class T>
  T prefixSumSerial(T *data, intT s, intT e)
  {
    T res = 0;
    for (intT i = s; i < e; ++i)
    {
      res += data[i];
      data[i] = res - data[i];
    }
    return res;
  }

  template <class T>
  void addSerial(T *data, intT s, intT e, T val)
  {
    for (intT i = s; i < e; ++i)
      data[i] += val;
  }

  /**
   * 对data的元素求prefixsum
   * data[0]=0
   * data[i]=data[0]+...+data[i-1]
   * T:return 所有数的和
   */
  template <class T>
  T prefixSum(T *data, intT s, intT e)
  {
    intT l = nblocks(e - s, _SCAN_BSIZE);
    if (l <= 1)
      return prefixSumSerial(data, s, e);
    T *sums = newA(T, l);
    blocked_for(i, s, e, _SCAN_BSIZE,
                sums[i] = prefixSumSerial<T>(data, s, e););
    T res = prefixSumSerial(sums, 0, l);
    blocked_for(i, s, e, _SCAN_BSIZE,
                addSerial(data, s, e, sums[i]););
    return res;
  }

  /** 
   * 将f中的元素中从s到t的元素按照F1布尔函数打包。对于i\in [s,t)，f(i)返回第i个元素
   * 即i从[s,t)，如果F1[i]=True,将f(i)元素顺次放入Out
   * 这个pack函数是最基本的，其他几个重载函数有更简洁的接口
   */
  template <class ET, class intT, class F>
  _seq<ET> pack(ET *Out, bool *Fl, intT s, intT e, F f)
  {
    intT l = nblocks(e - s, _F_BSIZE);
    if (l <= 1)
      return packSerial(Out, Fl, s, e, f);
    intT *Sums = newA(intT, l);
    //sums[i]放入F1第i个block的非0元素
    blocked_for(i, s, e, _F_BSIZE, Sums[i] = sumFlagsSerial(Fl + s, e - s););
    intT m = plusScan(Sums, Sums, l);
    if (Out == NULL)
      Out = newA(ET, m);
    //以block为单位进行pack（而不是以单个元素为单位pack）
    blocked_for(i, s, e, _F_BSIZE, packSerial(Out + Sums[i], Fl, s, e, f););
    free(Sums);
    return _seq<ET>(Out, m);
  }

  template <class ET, class intT, class F>
  void splitSerial(ET *OutFalse, ET *OutTrue, bool *Fl, intT s, intT e, F f)
  {
    intT kT = 0;
    intT kF = 0;
    for (intT i = s; i < e; i++)
      if (Fl[i])
        OutTrue[kT++] = f(i);
      else
        OutFalse[kF++] = f(i);
  }

  // Given a boolean array, splits so false (0) elements are at the bottom
  // and true (1) elements are at the top of the output (of lenght e-s).
  // As usual s is a start index, e is an end index and
  // f is a function of type [s,e-1) -> ET
  template <class ET, class intT, class F>
  int split(ET *Out, bool *Fl, intT s, intT e, F f)
  {
    intT l = nblocks(e - s, _F_BSIZE);
    intT *SumsTrue = newA(intT, l);
    blocked_for(i, s, e, _F_BSIZE, SumsTrue[i] = sumFlagsSerial(Fl + s, e - s););
    intT numTrue = plusScan(SumsTrue, SumsTrue, l);
    intT numFalse = (e - s) - numTrue;
    ET *OutTrue = Out + numFalse;
    blocked_for(i, s, e, _F_BSIZE,
                splitSerial(Out + _F_BSIZE * i - SumsTrue[i],
                            OutTrue + SumsTrue[i],
                            Fl, s, e, f););
    free(SumsTrue);
    return numFalse;
  }

  /**
  * 将f中的元素中从s到t的元素按照F1布尔函数打包。对于i\in [s,t)，f(i)返回第i个元素.
  * 打包后返回两个顺序数据对象pair<seq1,seq2>,
  * seq1打包了所有F1[i]=True的元素
  * seq2打包了所有F2[i]=True的元素
  */
  template <class ET, class intT, class F>
  pair<_seq<ET>, _seq<ET>> pack2(ET *Out, bool *Fl1, bool *Fl2,
                                 intT s, intT e, F f)
  {
    intT l = nblocks(e - s, _F_BSIZE);
    intT *Sums1 = newA(intT, l);
    intT *Sums2 = newA(intT, l);
    blocked_for(i, s, e, _F_BSIZE,
                Sums1[i] = sumFlagsSerial(Fl1 + s, e - s);
                Sums2[i] = sumFlagsSerial(Fl2 + s, e - s););
    intT m1 = plusScan(Sums1, Sums1, l);
    intT m2 = plusScan(Sums2, Sums2, l);
    ET *Out1;
    ET *Out2;
    if (Out == NULL)
    {
      Out1 = newA(ET, m1);
      Out2 = newA(ET, m2);
    }
    else
    {
      Out1 = Out;
      Out2 = Out + m1;
    }
    blocked_for(i, s, e, _F_BSIZE,
                packSerial(Out1 + Sums1[i], Fl1, s, e, f);
                packSerial(Out2 + Sums2[i], Fl2, s, e, f););
    free(Sums1);
    free(Sums2);
    return pair<_seq<ET>, _seq<ET>>(_seq<ET>(Out1, m1), _seq<ET>(Out2, m2));
  }

  /**
   * 该pack函数区别于上述函数在于用bit-vector，即F1来表示两种状态
   * 如果F1[i]=0,将f(i)打包在seq_的第一项
   * 如果F1[i]=1,将f(i)打包在seq_的第二项
   * 注意Out中是pair<_seq1,seq2>的容器
   */
  // Custom pack2 to be used with bitvector as flags (used for example for wavelet trees)
  template <class ET, class intT, class F>
  pair<_seq<ET>, _seq<ET>> pack2(ET *Out, long *Fl, intT s, intT e, F f)
  {
    // If interval empty
    if (s >= e)
      return pair<_seq<ET>, _seq<ET>>(_seq<ET>(Out, 0), _seq<ET>(Out, 0));
    intT l = nblocks(e - s, _F_BSIZE);
    intT *Sums1 = newA(intT, l);
    intT *Sums2 = newA(intT, l);
    blocked_for(i, s, e, _F_BSIZE,
                Sums2[i] = sumBitFlagsSerial(Fl, s, e); // count ones
                Sums1[i] = (e - s - Sums2[i]););        // calculate zeros
    intT m1 = plusScan(Sums1, Sums1, l);
    intT m2 = plusScan(Sums2, Sums2, l);
    ET *Out1;
    ET *Out2;
    if (Out == NULL)
    {
      Out1 = newA(ET, m1);
      Out2 = newA(ET, m2);
    }
    else
    {
      Out1 = Out;
      Out2 = Out + m1;
    }
    blocked_for(i, s, e, _F_BSIZE,
                packSerial01(Out1 + Sums1[i], Out2 + Sums2[i], Fl, s, e, f););
    //packSerial0(Out1+Sums1[i], Fl, s, e, f);
    //packSerial1(Out2+Sums2[i], Fl, s, e, f););
    free(Sums1);
    free(Sums2);
    return pair<_seq<ET>, _seq<ET>>(_seq<ET>(Out1, m1), _seq<ET>(Out2, m2));
  }

  template <class ET, class intT>
  intT pack(ET *In, ET *Out, bool *Fl, intT n)
  {
    return pack(Out, Fl, (intT)0, n, getA<ET, intT>(In)).n;
  }

  template <class ET, class intT>
  intT split(ET *In, ET *Out, bool *Fl, intT n)
  {
    return split(Out, Fl, (intT)0, n, getA<ET, intT>(In));
  }

  template <class ET, class intT>
  pair<intT, intT> pack2(ET *In, ET *Out, bool *Fl1, bool *Fl2, intT n)
  {
    pair<_seq<ET>, _seq<ET>> r;
    r = pack2(Out, Fl1, Fl2, (intT)0, n, getA<ET, intT>(In));
    return pair<intT, intT>(r.first.n, r.second.n);
  }

  // Custom pack which takes an input and one flag array and
  // puts all elements where 0 is set on the left side and
  // all the other elements on to the right
  template <class ET, class intT>
  intT pack2Bit(ET *In, ET *Out, long *Flags, intT s, intT e)
  {
    pair<_seq<ET>, _seq<ET>> r;
    r = pack2(Out, Flags, s, e, getA<ET, intT>(In));
    return r.first.n;
  }

  //打包In的n个元素到seq中
  template <class ET, class intT>
  _seq<ET> pack(ET *In, bool *Fl, intT n)
  {
    return pack((ET *)NULL, Fl, (intT)0, n, getA<ET, intT>(In));
  }

  //打包[0,n)的index到Out中，返回Out长度
  //即In[i]=i，bool[i]指示index i是否被打包
  template <class intT>
  intT packIndex(intT *Out, bool *Fl, intT n)
  {
    return pack(Out, Fl, (intT)0, n, utils::identityF<intT>()).n;
  }

  //打包[0,n)的index到seq线性类型中返回
  template <class intT>
  _seq<intT> packIndex(bool *Fl, intT n)
  {
    return pack((intT *)NULL, Fl, (intT)0, n, utils::identityF<intT>());
  }

  template <class ET, class intT, class PRED>
  intT filterSerial(ET *In, ET *Out, intT n, PRED p)
  {
    intT k = 0;
    for (intT i = 0; i < n; i++)
      if (p(In[i]))
        Out[k++] = In[i];
    return k;
  }

  /**
 * In中满足函数p的元素过滤到Out中，函数返回过滤后Out的长度
 * p:ET -> {0,1}
 */
  template <class ET, class intT, class PRED>
  intT filter(ET *In, ET *Out, intT n, PRED p)
  {
    if (n < _F_BSIZE)
      return filterSerial(In, Out, n, p);
    bool *Fl = newA(bool, n);
    parallel_for(intT i = 0; i < n; i++) Fl[i] = (bool)p(In[i]);
    intT m = pack(In, Out, Fl, n);
    free(Fl);
    return m;
  }

  /**
   *  In中满足函数p的元素过滤到Out中，函数返回过滤后Out的长度
   * F1 是一个长度为n的已经分配好的bool数组
   */
  // Avoids reallocating the bool array
  template <class ET, class intT, class PRED>
  intT filter(ET *In, ET *Out, bool *Fl, intT n, PRED p)
  {
    if (n < _F_BSIZE)
      return filterSerial(In, Out, n, p);
    parallel_for(intT i = 0; i < n; i++) Fl[i] = (bool)p(In[i]);
    intT m = pack(In, Out, Fl, n);
    return m;
  }

  /**
   *   In中满足函数p的元素过滤出来，放入_seq中放回
   *  实际就是filter and pack
   */
  template <class ET, class intT, class PRED>
  _seq<ET> filter(ET *In, intT n, PRED p)
  {
    bool *Fl = newA(bool, n);
    parallel_for(intT i = 0; i < n; i++) Fl[i] = (bool)p(In[i]);
    _seq<ET> R = pack(In, Fl, n);
    free(Fl);
    return R;
  }

  //如果结果很少（In中很多都元素都被淘汰），该filter算法更快
  // Faster for a small number in output (about 40% or less)
  // Destroys the input.   Does not need a bool array.
  template <class ET, class intT, class PRED>
  intT filterf(ET *In, ET *Out, intT n, PRED p)
  {
    intT b = _F_BSIZE;
    if (n < b)
      return filterSerial(In, Out, n, p);
    intT l = nblocks(n, b);
    b = nblocks(n, l);
    intT *Sums = newA(intT, l + 1);
    {
      parallel_for(intT i = 0; i < l; i++)
      {
        intT s = i * b;
        intT e = min(s + b, n);
        intT k = s;
        for (intT j = s; j < e; j++)
          if (p(In[j]))
            In[k++] = In[j];
        Sums[i] = k - s;
      }
    }
    intT m = plusScan(Sums, Sums, l);
    Sums[l] = m;
    {
      parallel_for(intT i = 0; i < l; i++)
      {
        ET *I = In + i * b;
        ET *O = Out + Sums[i];
        for (intT j = 0; j < Sums[i + 1] - Sums[i]; j++)
        {
          O[j] = I[j];
        }
      }
    }
    free(Sums);
    return m;
  }

} // namespace sequence
#endif // _A_SEQUENCE_INCLUDED
